…
To reproduce this workflow, first clone the repository to your local machine or cluster:
# git clone https://github.com/mandymejia/BayesianBrainMapping-Templates.git
# cd BayesianBrainMapping-TemplatesThis section initializes the environment by loading required packages, setting analysis parameters, and defining directory paths.
Important: Before running the workflow, you must
update the following variables in 0_setup.R to match your
local or cluster environment:
dir_project, dir_data,
dir_slate, dir_personal
HCP_restricted_fname (path to the restricted HCP CSV
if you have access to it)
wb_path (location of the CIFTI Workbench on your
system)
github_repo_dir <- getwd()
src_dir <- file.path(github_repo_dir, "src")
source(file.path(src_dir, "0_setup.R"))## Using this Workbench path: '/Users/nohelia/Downloads/workbench/bin_macosxub/wb_command'.Before estimating group-level templates, we apply a set of filtering steps to select a high-quality subject sample. These steps ensure that the final templates are based on reliable, representative data. The filtering steps include:
We begin by filtering subjects based on the duration of valid fMRI
data after motion correction. For each subject, and for each session
(REST1, REST2) and encoding direction
(LR, RL), we compute Framewise Displacement
(FD) using the fMRIscrub package. FD is calculated from the
Movement_Regressors.txt file available in the HCP data for
each subject, encoding and session.
A volume is considered valid if it passes an FD threshold, and a subject is retained only if both sessions in both encodings have at least 10 minutes (600 seconds) of valid data.
The final subject list includes only those who passed the filtering
criteria in both LR and RL encodings. This list is referred to as the
combined list and is the one used throughout this
project.
# source(file.path(src_dir,"1_fd_time_filtering.R"))During this step, an FD summary table is generated with the following columns:
subject: HCP subject ID
session: REST1 or REST2
encoding: LR or RL
mean_fd: mean framewise displacement
valid_time_sec: total duration of valid data in seconds
# Read FD summary
fd_summary <- read.csv("~/Documents/StatMIND/Data/fd_summary.csv")
# Display the first 4 rows
knitr::kable(head(fd_summary, 4), caption = "First rows of FD summary table")| X | subject | session | encoding | mean_fd | valid_time_sec |
|---|---|---|---|---|---|
| 1 | 100206 | REST1 | LR | 0.1017240 | 858.24 |
| 2 | 100206 | REST2 | LR | 0.1361220 | 858.96 |
| 3 | 100206 | REST1 | RL | 0.0698779 | 864.00 |
| 4 | 100206 | REST2 | RL | 0.0824894 | 863.28 |
As shown above, subject 100206 qualifies for further analysis because each of the four sessions (REST1/REST2 × LR/RL) contains at least 600 seconds of valid data.
The script is currently designed to filter based on valid time only, but it can be easily adapted to apply additional constraints such as maximum mean FD thresholds if desired (e.g., mean_fd < 0.2).
In the final step of the subject filtering pipeline, we balance sex across age groups to reduce potential demographic bias in template estimation.
For the combined list of valid and unrelated subjects,
we:
Subset the HCP unrestricted demographics to include only those subjects.
Split subjects by age group and examine the sex distribution within each group.
If both sexes are present but imbalanced, we randomly remove subjects from the overrepresented group to achieve balance.
Note: If you are not applying the unrelated subject filtering step
(3.2), you can modify the code to subset based on
valid_combined_subjects_FD instead of
valid_combined_subjects_unrelated.
The final list of valid subjects is saved as:
valid_combined_subjects_balanced.csv
valid_combined_subjects_balanced.rds (used in the
template estimation step)
#source(file.path(src_dir,"3_balance_age_sex.R"))In this step, we load and preprocess a group-level cortical
parcellation to be used to estimate template in the next step.
Specifically, we use the Yeo 17-network parcellation
(Yeo_17) and perform the following operations:
Simplify the labels by collapsing hemisphere-specific naming and removing subnetwork identifiers, grouping regions by their main network.
Create a new dlabel object that maps each vertex to
its corresponding network.
Mask out the medial wall to exclude non-cortical regions from analysis.
The resulting parcellation is saved as
Yeo17_simplified_mwall.rds.
# source(file.path(src_dir,"4_parcellations.R"))We can visualize the Yeo17 networks and their corresponding labels:
# Load libraries
library(ciftiTools)
library(rgl)
rgl::setupKnitr()
# Load the parcellation
yeo17 <- readRDS(file.path(dir_data, "Yeo17_simplified_mwall.rds"))
yeo17 <- add_surf(yeo17)
view_xifti_surface(
xifti = yeo17,
widget = TRUE,
title = "Yeo17 Network Parcellation",
legend_ncol = 6,
legend_fname = "yeo17_legend.png",
)